Drag force scaling for penetration into granular media 
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Impact dynamics is measured for spherical and cylindrical projectiles of many different densities 
dropped onto a variety non-cohesive granular media. The results are analyzed in terms of the 
material-dependent scaling of the inertial and frictional drag contributions to the total stopping 
force. The inertial drag force scales similar to that in fluids, except that it depends on the internal 
friction coefficient. The frictional drag force scales as the square-root of the density of granular 
medium and projectile, and hence cannot be explained by the combination of granular hydrostatic 
pressure and Coulomb friction law. The combined results provide an explanation for the previously- 
observed penetration depth scaling. 

PACS numbers: 45.70.Ht, 45.70.Cc, 83.80.Fe, 89.75.Da 



Collision is one of the most fundamental processes in 
nature, and can be exploited to uncover the basic physics 
of systems ranging from planets to elementary particles. 
Impact of projectiles into a pack of grains has been of in- 
creasing interest for highlighting and elucidating the un- 
usual mechanical properties of granular materials. This 
includes studies of crater morphology [IH3], penetration 
depth [11], dynamics HHIEj, boundary effects [TBI [T^ 
[23] and packing- fraction effects [HI |23l [24] . One finding 
is that the penetration depth scales as d ^ Dp^'^H^^^ 
where Dp is projectile diameter and H is the total drop 
distance [4j |5]- Granular impact is also important for 
military and industrial applications |25[ 126] . and cone 
penetration tests are used for in-situ soil characteriza- 
tion [13 [Ml- 

In a previous study [13] we measured the impact dy- 
namics of a 2.54 cm steel ball onto a packing of glass 
beads. The stopping force was found to be the sum 
of an inertial drag, proportional to the square of the 
speed V, and a frictional drag, proportional to depth z. 
This has since been supported by several other experi- 
ments [TU [TBtiTSl [2Tti23] [29] . The equation of motion 
during impact is thus 



ma = —mg + mv /di + k\z\, 



(1) 



where m is projectile mass, g — 980 cm/s^, and di 
and k are materials parameters expressed with units of 
a length and a spring constant, respectively. By using 
ma = dK/dz and an integrating factor, as in Ref. [TU] . 
this can be solved for speed versus depth: 
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where vq is the impact speed at z = 0. The final pene- 
tration depth d is then given by the limit of v — > as 
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where W{x) is the Lambert W- function. An additional 
constant stopping force /o [SldH] can be included in these 
expressions by replacing g with g — fo/m. 

Some open questions are how di and k scale with the 
materials properties of the projectile and the granular 
packing, and how this conspires to give d ^ Dp H^/^. 
For inertial drag, exactly as for hydrodynamic drag at 
high Reynolds number Re, momentum transfer gives 
the expectation mv'^/di ~ Apgv'^ where A is the pro- 
jected projectile area and Pg is the mass density of the 
granular medium. For frictional drag, the combina- 
tion of hydrostatic pressure and Coulomb friction gives 
k\z\ ~ figpgA\z\ where the internal frictional coefficient 
is /i = tan^r and 9r is the angle of repose. The scaling 
of di and k would thus be 



di/Dp - Pp/Pg, 
kDp/mg - pPg/Pp, 



(4) 
(5) 



(3) 



where Dp and pp correspond to diameter and density 
of projectile, respectively. Here, we test the speed versus 
depth prediction of Eq. (pi) , and we compare fitted values 
of di and k with the above expectations. While neither 
turns out to be quite correct, we connect the results to 
the observed penetration depth scaling. 

Our basic experimental setup is identical with previous 
work [13]. The velocity of the projectile v{t) at time t is 
computed by particle-image velocimetry applied to fine 
stripes on a vertical rod glued to the top of the projectile. 
The system has 20 /iS temporal resolution and 100 nm 
spatial resolution, which is fine enough to compute posi- 
tion and acceleration from v(t). The primary difference 
from Ref. [13] is that now we vary both the projectiles and 
the granular medium. We begin with 0.35 mm diameter 
glass beads, prepared to a reproducible random packing 
state by slowly turning off a fluidizing up-flow of air. Into 
this we drop a wide variety of projectiles, listed in Ta- 
ble |lj from free-fall heights between to 85 cm. The first 
type is spheres. The second type is aluminum cylinders, 
which are dropped lengthwise with the axis horizontal 



Projectile 




Pv (g/cm^*) 


Dp (cm) 


Lp (cm) 


Hollow PP ball 


0.51 


2.54 


- 


Wood ball 




0.95 


2.54 


- 


Delrin ball 




1.65 


2.54 


- 


PTFE ball 




2.46 


2.54 


- 


Steel ball 




7.96 - 159 


0.3175 - 5.08 


- 


Tungsten carbide ball 


15.3 


2.54 


- 


Aluminum 


cylinder 


2.89 - 4.26 


0.635 -1.27 


5.08 - 15.24 



TABLE I: Projectile properties. The steel sphere diameters 
are Dp = 5.08, 2.54, 1.27, 0.635, and 0.3175 cm. The alu- 
minum cylinders diameters are Dp = 1.27 and 0.635 cm; each 
have lengths Lp = 5.08, 10.16, and 15.24 cm. The density pp 
is projectile plus rod mass divided by projectile volume. 



and parallel to the surface of the granular medium. In 
both cases the effective density is given by projectile plus 
rod mass divided by projectile volume. The length of the 
cylinders is varied, but we find that it does not affect the 
dynamics or final penetration depths. 

Example speed versus position data are plotted in 
Fig.fllfor Dp — 2.54 cm diameter projectiles of four differ- 
ent densities, dropped onto the glass bead packing, with 
initial impact speeds ranging from zero to 400 cm/s. For 
slow impacts, the speed first increases and then decreases 
with depth. For faster impacts, the speed versus depth 
curves gradually change from concave down to concave 
up. Generally, there is a rapid decrease of speed to zero at 
the final penetration depth. We obtain good fits to these 
complex trajectories by adjusting k and di, as shown in 
Fig. [T] The displayed level of agreement is typical for all 
projectiles, including the cylinders. 

Unfortunately a good simultaneous fit for a given pro- 
jectile to a single pair of k and di values at all impact 
speeds can be obtained only for denser projectiles. For 
the less dense projectiles, the fitting parameters become 
constant only at high impact speeds. Then the same 
values for k and di are obtained as from the analysis 
method of Ref. [B]. This holds roughly for d > Dj,/2 
and Pp > 2pg as shown in Fig. 2] We speculate that for 
low density projectiles and small impact speeds, the pen- 
etration can be shallow enough that the detailed shape 
of projectile must be taken into account |30| . In this 
regime surface flows and surface roughness could plays 
a role, though identical penetration behavior was found 
earlier for slick and tacky projectiles of the same size and 
density [4j. We note that including a constant force /o 
as a third free parameter does not noticeably change the 
fits; the largest fit value is fo/m — 2.4 cm/s^, which is 
small compared to g. We note, too, that /o cannot be 
chosen such that k and di become constant. For the rest 
of the paper, we restrict attention to conditions where the 
deduced k and di values do not depend on impact speed 
and hence can be considered as materials parameters. 

We now compare the fitting parameters with Eqs (4]l5 1. 
The expectation for di is based on momentum transfer. 
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FIG. 1: (color online) Example velocity versus depth data 
for Dp — 2.54 cm diameter spheres. The impact begins at 
z — 0, and proceeds downward in the —z direction. The 
black dashed curves are a simultaneous fit to Eq. (pi, where 
a single pair of k and di values is found for different initial 
impact speeds. The grey dotted curves are also fits to Eq. pj, 
but where k and di are adjusted for each impact speed and 
hence do not represent well-defined materials parameters. 



so that the inertial drag force is mv'^/di ^ ApgV^ just like 
an object moving in a fluid at high Reynolds number. For 
spheres and cylinders, the characteristic length is thus 
di ~ m/Apg and can be written as di ^ D'pPp/ pg if we 
take Dp to be 1 times diameter for spheres and Svr/S 
times diameter for cylinders (the cylinder length drops 
out of the ratio m/ A). For a unified analysis, we therefore 
divide the fitted value of di by D' and plot versus Pp in 
Fig. [3^. We find that all data, including the cylinder 
data, collapse nicely onto a power-law with the expected 




FIG. 2: (color online) Scatter plot of penetration depth versus 
projectile density, scaled respectively by projectile diameter 
and grain density. An open circle is plotted for conditions 
such that the k and di fitting parameters are constant; a 
cross is plotted otherwise. Blue is for glass beads; green is for 
rice; pink is for beach sand; yellow is for sugar. 



slope of one; i.e. d/D' ^ p^ holds as per expectation. 

For the fitting parameter k that sets the quasi-static 
friction force, k\z\, we now make a similar compari- 
son with expectation by plotting kD' /mg versus pp in 
Fig. ^p. All the data, including the cylinder, again col- 
lapse nicely to a power law in projectile density. However, 
the expected power-law kDp/mg ^ ^/ Pp is clearly wrong. 
The data are instead consistent with kDp/mg ~ 1/ ,Jpp. 
Therefore, the nature of the quasistatic frictional drag is 
different from the simple combination of Coulomb fric- 
tion and hydrostatic pressure. 

To investigate this further, we perform a second series 
of experiments where £)p=2.54 cm diameter steel spheres 
are dropped into rice, beach sand, and sugar (with ma- 
terials properties listed in Table |ll]). As in Fig. [I] speed 
vs position data for a range of impact speeds are fit to 
Eq. ([2]) to obtain values of di and k. Based on Fig. |3J 
and assuming that the x-axis of this figure is correctly 
made dimensionless by dividing out the bulk density pg 
of the granular medium, the observed scalings so far are 
di/D'j, ~ {pp/pg) and kD'^/mg ~ {Ppl Pg)'^'"^ ■ There- 
fore, we divide out this density dependence and plot it in 
Fig. |4] versus the only remaining material property of the 
grains: the internal friction coefficient p = tan^,. where 
Or is the draining angle of repose. In this figure the data 
for glass beads, from Fig.|3J all lie at /x = 0.45. The range 
of p values is less than a factor of two, but to within 
uncertainty the scaled di and k parameters collapse to 
power-laws in p. For the quasistatic frictional drag co- 
efficient we find k ^ p, which is expected for Coulomb 
friction. For the speed-squared inertial drag coefficient, 
we find di ^ l//i. Therefore the inertial drag force is 
proportional to /i, which could correspond to an added- 
mass effect whereby the volume of grains boosted to the 
projectile speed grows in proportion to p. This is unlike 
the case of simple fluids, where the fluid flow and the 
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FIG. 3: (color online) (a) Dimensionless inertial drag length 
scale di/D'p, and (b) dimensionless quastistatic drag coeffi- 
cient kD'p/mg, versus projectile density, for impact into glass 
beads. Here Dp is the effective projectile diameter; m is its 
total mass, including rod; pp is m divided by projectile vol- 
ume; g — 980 cm/s^. Open symbols are for spheres made of 
hollow PP (>), wood (<), PTFE (v), delrin (A), steel (Q), 
and tungsten carbide (D); closed diamonds are for aluminum 
cylinders. Symbol sizes increase with Dp, with values given 
in Table [l] The solid gray lines denote power-laws as labeled. 
The dashed line in (b) is the expected scaling, kD'p/mg ~ Pp^. 



Granular Material 


Pa 


(g/cm^) 


Or 


Grain size (mm) 


Glass beads 




1.45 


24° 


0.25 - 0.35 


Rice 




0.77 


35° 


2x8 


Beach sand 




1.51 


36° 


0.2 - 0.8 


Sugar 




0.89 


40° 


0.4 - 0.7 



TABLE II: Granular media properties: pg is bulk density; 0,- 
is drainage angle of repose; size is the range of grain diameters, 
except for rice where it is the length of short and long axes. 



inertial drag force at high Re depend only on the density 
of the fluid. 

As an alternative analysis for the p dependence of di , 



one could imagine that an inertial drag force 



'A 



loads the contacts and gives an additional friction force 
of p times this loading. Then the total velocity-squared 
force is mv^ /di ^ {l+ap)pgv'^A, which gives di ex 1/(1 + 
ap). This reasonably flts the data, as shown in Fig. [4k 
with a = 2.2 ± 0.6. The residuals are smaller for the 
power-law form. 

Combining the power-laws in Figs. |3]|4] and the ac- 
tual numerical prefactors, the inertial and frictional drag 



Q_ 

" o. 

D 



0.1 



£ 10 



60 

S 



_ 




1 1 1 1 1 


_ 


^ -1 


T 






: - - _ _ 


^ 


h 


^0 


1 

/ 


1 1 1 1 1 


1 

1 1 


, , , , 1 


^ 


1 
- (b) 


1 1 


1 1 1 1 1 


- 






o^- 


^^ 


1 


- 


^ 


^^o 






1 
1 1 


1 1 1 1 1 


- 



10 



FIG. 4: (color online). (a) {di/D'j,){pg/ pp) and (b) 
{kD'p/mg){pp/ pgY'"^ , versus p = tanS,., where 9r is the re- 
pose angle of the granular medium. The symbols &t p = 0.45 
are for glass beads, with the same symbol codes as in Fig. [S] 
The colored circles at increasing p are for rice (green), beach 
sand (pink), and sugar (yellow). The gray lines indicate 
power-law behavior. The dashed curve in (a) is a fit to 
di oc 1/(1 + ap), giving a = 2.2 ± 0.6. 



coefficients are altogether found by measurements for a 
range of projectiles and grains to be consistent with 

(6) 

(7) 



di/i?; = (0.25/m)(Pp/p<,), 
kD'/mg = 12fi{pg/ pp)^^"^ . 



The two main differences from the simple expectation of 
Eqs. (|4]l5|, are the factor of 1/p in di and the density 
ratio exponent of 1/2 rather than 1 in k. These results 
may be inspected differently by re-writing the equation 
of motion as 

ma — —mg + 2.7 ppgV^ A + 8.0p {ppPg) g\z\A. (8) 

Note that the numerical prefactor for the depth- 
dependent frictional drag is significantly larger than 
unity, and that y/PpPg is larger than pg for dense pro- 
jectiles. For both reasons, frictional drag exceeds the 
value expected from hydrostatic pressure and Coulomb 
friction. One might have expected the opposite effect, 
either by a decrease in contact area between projectile 
and grains due to ejection of grains or by fluidization of 
the grains from the motion of the projectile. Our re- 
sults instead appear to indicate that frictional contacts 
are loaded by the motion of the projectile, so that the 
medium is stronger than set by gravity alone. Such be- 
havior is not seen for the fast horizontal rotation of bars 
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FIG. 5: (color online). Final penetration depth scaling plot. 
All trials data are shown and they reasobaly agree with the 
empirical scahng d = 0.14/x"^(pp/pg)^/^Dp^/^//^/^ (dotted 
line) [H |5] . Colored curves show the force law predictions. 
They are also close to data and the empirical scaling. 



P5] , where the depth- and speed-dependent forces are 
easily disentangled, and warrants further attention. 

As a final test we now compare data for the final pen- 
etration depth d with Eq. (Is]) . Prior observations [4j [5] 
are consistent with the empirical form 



d = 0.14/i-i(p,/p,)V2i^'V3^i/3. 



(9) 



Thus in Fig. [s] we plot penetration depth data for all 
trials versus the quantity {pp/ pgY^'^Dp H^/^ from the 
right-hand side of this expression. This collapses our new 
data to within the experimental uncertainty, including 
that for the cylinder, to the line y = 0.14a; and hence 
shows agreement with prior work. However Eqs. ( 3|8 ) 
do not predict perfect power-law behavior. To compare 
with data, we first calculate the geometric mean of each 
of the five variables H, Dp, pp, pg, and p, over the range 
of experimental conditions employed here. The predicted 
penetration depth for these mean values is shown by a 
single red open circle in Fig. [51 Then we vary each of vari- 
ables, one at a time, keeping all others fixed at their mean 
value. The resulting five predicted penetration depth 
curves are included in Fig. [5] They are not identical, 
or even perfect power laws, but are all close together and 
in fair agreement with the data. Thus the empirical pen- 
etration depth data scaling is satisfactorily understood 
in terms of the nature of the stopping force exerted by 
the medium onto the projectile. 

In conclusion, we developed an exact solution of Eq. (fTl) 
for the dynamics of penetration, and we conducted a 
wide range of granular impact experiments to elucidate 
the materials dependence of the inertial and frictional 
contributions to the stopping force. The final equation 
of motion, Eq. (Is]), is roughly consistent with the em- 
pirical scaling of penetration depth versus drop height 
and materials parameters. However it is not consistent 
with the apparently simplistic expectation of Eqs. (|4]l5|. 
So there is unanticipated physics, yet to be understood, 



which could possibly arise from motion-loading of fric- 
tional contacts or from granular flow fields that depend 
on the internal friction coefficient. 

This work was supported by the JSPS Postdoctoral 
Fellowships for Research Abroad (HK) and the National 
Science Foundation through Grant Nos. DMR-0704147 
and DMR-1305199. 



[1 
[2 
[3 

[4; 

[5 
[6 

[7: 

[10 

[11 
[12 

[13 



J. C. Amato and R. E. Williams, Am. J. Phys. 66, 141 

(1998). 

A. M. Walsh, K. E. HoUoway, P. Habdas, and J. R. 

de Bruyn, Phys. Rev. Lett. 91, 104301 (2003). 

X. J. Zheng, Z. T. Wang, and Z. G. Qiu, Eur. Phys. J. E 

13, 321 (2004). 

J. S. Uehara, M. A. Ambroso, R. P. Ojha, and D. J. 

Durian, Phys. Rev. Lett. 90, 194301 (2003). 

M. A. Ambroso, C. R. Santore, A. R. Abate, and D. J. 

Durian, Phys. Rev. E 71, 051305 (2005). 

J. R. de Rruyn and A. M. Walsh, Can. J. Phys. 82, 439 

(2004). 

K. E. Daniels, J. E. Coppock, and R. P. Behringer, Chaos 

14, S4 (2004). 

M. P. Ciamarra, A. H. Lara, A. T. Lee, D. L Goldman, 

L Vishik, and H. L. Swinney Phys. Rev. Lett. 92, 194301 

(2004). 

D. Lohse, R. Rauhe, R. Bergmann, and D. van der Meer, 

Nature 432, 689 (2004). 

M. A. Ambroso, R. D. Kamien, and D. J. Durian, Phys. 

Rev. E 72, 041305 (2005). 

M. Hou, Z. Peng, R. Liu, K. Lu, and C. K. Chan, Phys. 

Rev. E 72, 062301 (2005). 

L. S. Tsimring and D. Volfson, in Powders and Grains 

2005, edited by R. Garcia- Rojo, H. Herrmann, and S. Mc- 

Namara (Balkema, Rotterdam, 2005), p. 1215. 

H. Katsuragi and D. J. Durian, Nature Phys. 3, 420 



[14 
[15 
[16 

[17 
[18 
[19 
[20 
[21 
[22 

[23' 
[24 
[25 
[26 
[27 
[28 
[29 
[30 



(2007). 

D. L Goldman and P. B. Umbanhowar, Phys. Rev. E 77, 
021308 (2008). 

P. B. Umbanhowar and D. I. Goldman, Phys. Rev. E 82, 

010301 (2010). 

F. Pachoco- Vazquez, G. A. Caballero-Robledo, J. M. 

Solano-Altamirano, E. Altshuler, A. J. Batista-Leyva, 

and J. C. Ruiz-Suarez, Phys. Rev. Lett. 106, 218001 

(2011). 

A. H. Clark, L. Kondic, and R. P. Behringer, Phys. Rev. 

Lett. 109, 238302 (2012). 

A. H. Clark and R. P. Behringer, Europhys. Lett. 101, 

64001 (2013). 

J. F. Boudet, Y. Amarouchene, and H. Kellay, Phys. Rev. 

Lett. 96, 158001 (2006). 

E. L. Nelson, H. Katsuragi, P. Mayor, and D. J. Durian, 
Phys. Rev. Lett. 101, 068001 (2008). 

A. Seguin, Y. Bertho, and P. Gondret, Phys. Rev. E 78, 

010301 (2008). 

S. von Kann, S. Joubaud, G. A. Caballero-Robledo, 

D. Lohse, and D. van der Meer, Phys. Rev. E 81, 041306 

(2010). 

H. Torres, A. Gonzalez, G. Sanchcz-Colina, J. C. Drake, 

and E. Ahshuler, Rev. Cub. Fis. 29, 1E45 (2012). 

J. R. Royer, B. Conyers, E. L Corwin, P. J. Eng, and 

H. M. Jaeger, Europhys. Lett. 93, 28008 (2011). 

D. J. Roddy, R. O. Pepin, and R. B. Merril, eds., Impact 

and explosion cratering (Pergamon Press, NY, 1978). 

J. A. Zukas, ed.. High velocity impact dynamics (Wiley, 

NY, 1990). 

H. S. Yu and J. K. MitcheU, J. Geotech. Eng. 124, 140 

(1998). 

M. Toiya, J. Hettinga, and W. Losert, Granular Matter 

9, 323 (2007). 

T. A. Brzinski and D. J. Durian, Soft Matter 6, 3038 

(2010). 

K. A. NewhaU and D. J. Durian, Phys. Rev. E 68, 

060301 (R) (2003). 



